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I Abstract 

We describe a molecular dynamics framework for the direct calculation of the short-ranged struc- 
tural forces underlying grain-boundary premelting and grain-coalescence in solidification. The 
method is applied in a comparative study of (i) a £9 (115) 120° twist and (ii) a £9 (110) {411} 



symmetric tilt boundary in a classical embedded-atom model of elemental Ni. Although both 



boundaries feature highly disordered structures near the melting point, the nature of the tempera- 

m 

ture dependence of the width of the disordered regions in these boundaries is qualitatively different. 
The former boundary displays behavior consistent with a logarithmically diverging premelted layer 
thickness as the melting temperature is approached from below, while the latter displays behavior 

o 

t-h featuring a finite grain-boundary width at the melting point. It is demonstrated that both types 

of behavior can be quantitatively described within a sharp-interface thermodynamic formalism 

>< 

involving a width- dependent interfacial free energy, referred to as the disjoining potential. The 
disjoining potential for boundary (i) is calculated to display a monotonic exponential dependence 
on width, while that of boundary (ii) features a weak attractive minimum. The results of this work 
are discussed in relation to recent simulation and theoretical studies of the thermodynamic forces 
underlying grain-boundary premelting. 
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INTRODUCTION 



At high homologous temperatures the atomic structure of a grain boundary often dis- 
plays pronounced disorder. In some cases this structural disorder can involve the formation 
of nanometer-scale intergranular films with liquid-like properties below the bulk melting 
point, a phenomenon commonly referred to as grain-boundary premelting. The interfacial 
thermodynamic driving forces underlying grain-boundary premelting are understood to be 
an important factor influencing grain coalescence behavior during the late stages of solidifi- 
cation (e.g., [lj). Specifically, when premelting is thermodynamically favored there exists an 
associated repulsive "disjoining pressure" which hinders the coalescence of two misoriented 
grains at nanometer-scale distances. In such cases a significant undercooling may be required 
for dendrite arms to merge to form solid "bridges" that are capable of sustaining thermal- 
contraction stresses without grain sliding or rupture. The quantitative characterization of 
grain-boundary disjoining forces is thus an important issue in the context of modeling the 
formation of solidification defects, known as "hot tears," which occur deep within the mushy 
zone during casting or welding pJH]. 

Despite its practical importance in the context of solidification defects, the magnitude and 
spatial extent of grain-boundary disjoining forces in metals and alloys remain incompletely 
understood. As in the case of surface premelting [5], these forces can in principle be probed 
experimentally through measurements of the extent of equilibrium premelting as a function 
of temperature near the bulk melting point. In comparison to surface premelting, however, 
the challenges inherent in characterizing the structure of "buried" internal interfaces at 
high homologous temperatures have significantly limited the number of direct experimental 
studies of grain-boundary premelting [6HT3]. This situation has provided motivation for 
several recent theoretical studies based on conventional phase-field [T4l - [TT] and phase-field- 
crystal (PFC) p~8j [19] methods, which have led to new insights into the rich variety of 
possible premelting behavior that may be exhibited by grain boundaries as a function of 
their bi-crystallography. 

Due to the difficulties inherent in validating theoretical models for grain-boundary pre- 
melting directly from experimental measurements, the authors recently proposed an in- 
dependent methodology for calculating grain-boundary disjoining forces from equilibrium 
molecular-dynamics (MD) simulations [20]. The technique represents an extension of the 
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numerous previous MD studies of grain-boundary premelting performed over the past three 
decades [2TH34] . The technique for calculating disjoining potentials by MD is based on an 
analysis of the equilibrium distribution of the widths of premelted intergranular films, which 
are related to the underlying disjoining forces through the fluctuation-dissipation theorem. 
The purpose of this paper is to describe the technical details surrounding the implementa- 
tion of this method, and to demonstrate its application in the study of two distinct classes 
of grain boundary premelting behavior. Specifically, we consider the premelting behavior of 
a high and intermediate energy boundary in elemental Ni where the widths of the premelted 
layers continuously increase or remain finite, respectively, as the melting point is approached 
from below. It is shown that the former behavior is consistent with an interfacial free energy 
that decreases exponentially with increasing width (w) of the premelted layer, while the 
latter behavior can be quantitatively modeled with a dependence of \& on w that features 
a weak attractive minimum at nanometer-scale widths. This minimum is accompanied by 
a grain boundary structure with alternating solid bridges and disordered regions that bears 
strong similarities to the type of structure found to be associated with such a minimum in 
a recent phase-field-crystal study [19]. 

The remainder of this paper is outlined as follows. In the next section we review a 
continuum, sharp-interface formalism for the thermodynamic properties underlying grain- 
boundary premelting and coalescence, based on the definition of the so-called "disjoining 
potential," the (negative) derivative of which is the disjoining pressure referred to above. 
Section III discusses the technical details underlying the calculation of the disjoining poten- 
tial by MD, as well as the details of the simulations undertaken in the current application 
of the method to the study of high-temperature grain boundaries in fee Ni. The results are 
presented in section IV, followed by a summary and discussion of the findings in the context 
of previous theoretical studies in section V. 

DISJOINING POTENTIAL 

The equilibrium width of a premelted grain-boundary is governed by a competition be- 
tween bulk and interfacial thermodynamic factors, which can be represented mathematically 
as follows: 

G{w) = AG f w + ^(w). (1) 
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In Eq.[TJ G(w) represents the total excess free energy of a premelted grain boundary of width 
w, which is composed of two contributions: AGf represents the bulk free energy difference 
between liquid and solid phases (positive below the melting point), per unit volume, while 
ty(w) represents the width-dependent interfacial free energy, which we refer to as the dis- 
joining potential. The function ty(w) takes the limits of jgb (the interfacial free energy of a 
hypothetical "dry" grain boundary) and 2^sl (twice the solid-liquid interfacial free energy) 
for small and infinite values of w, respectively. For intermediate values of the width ty(w) 
can display a complex dependence on w. 

The dependence of ty(w) on w is generally governed by distinct short and long-ranged 
contributions. An attractive interaction between solid-liquid interfaces (35J [36] arises due 
to dispersion forces which are dominant at large w, and are predicted to give rise to finite 
interfacial widths at Tm [35] . For systems where the wetting condition, jgb > ^Isl holds, a 
repulsive contribution to ty(w) arises from short-ranged structural interactions (\? sr ) associ- 
ated with the overlap of the density waves in the diffuse regions of the solid-liquid interfaces. 
Mean-field arguments [37J, as well as lattice-gas models (e.g., [28]). yield an exponentially 
decaying form for this short-ranged contribution: 

^ sr (w) = 2-fsL + A 7 exp[-w/5] (2) 

where A7 = ^qb — and 5 is an interaction length on the order of the atomic spacing. 
In the absence of long-ranged dispersion forces, insertion of Eq. [2] into Eq. [j] leads to the pre- 
diction of an equilibrium grain boundary width that diverges logarithmically as the melting 
temperature (Tm) is approached from below. In previous work [20] we have argued, based 
on previous estimates of the dispersion forces for surface premelting, that the structural con- 
tributions to the disjoining potential for grain-boundary premelting in metals are expected 
to dominate the long-ranged dispersion contributions for nanometer-scale grain-boundary 
widths. 

Recent theoretical results demonstrate that ^ sr (w) can generally display more complex 
dependencies on w than suggested by Eq. [2| Diffuse-interface phase field models [T1HT6] . 
which neglect dispersion forces and thus model \^ 5r directly, have shown that depending on 
grain-boundary misorientation and the detailed choice of interfacial thermodynamic param- 
eters, the premelting transition can exhibit either the continuous behavior associated with 
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Eq. [2j a hysteretic first-order character, or an intermediate behavior where the boundary 
width increases with increasing temperature but remains finite at T M (these three types of 
behavior are referred to as types 2, 1 and 3 in Ref. [34]). The PFC method was recently used 
to study grain boundary premelting for both three-dimensional body-centered-cubic (bcc) 
[T8] and two-dimensional hexagonal [19J systems. The latter study involved a systematic 
investigation of symmetric tilt boundaries as a function of misorientation. Well beyond a 
critical misorientation angle (where ^gb > ^Isl) ^ sr was found to exhibit a purely "repul- 
sive" behavior, monotonically decreasing with increasing w ) consistent with Eq.[2j However, 
well below the critical angle, ty sr exhibited an attractive minimum, giving rise to a finite 
area-averaged liquid-layer width at the melting temperature, reflecting the presence of local- 
ized premelting within the cores of the grain-boundary dislocations in low- angle boundaries. 

The possibility that the short-ranged contributions to the disjoining potential can exhibit 
an attractive minimum was considered in earlier theoretical studies of wetting transitions 
[38] . In these studies a double-exponential ansatz was used to model such behavior. For the 
purposes of the present study this double-exponential form can be written as follows: 

^ sr (w) = 2 1SL + [Ai exp (-w/Sr) - A 2 exp {-w/6 2 )], (3) 

where Ai and A 2 , which are both positive quantities, represent the strengths of repulsive 
and attractive exponential contributions with decay lengths Si and 5 2 , respectively. Al- 
though Eq. [3] has not been derived directly from a microscopic theory, it provides a useful 
phenomenological form to model the MD data, as shown below. This equation gives rise 
to a predicted value of the quantity Jgb — ^Isl equal to Ai — A 2 , where again jgb is the 
interfacial free energy of a hypothetical dry grain boundary given by the limit of w going to 
zero in Eq. [3| Choosing 5 2 > 5i, and if Ai > A 2 , as in the results given below, this disjoining 
potential is repulsive at short distances, attractive at large w, and features a minimum at a 
width w m given as: 

With the form for the disjoining potential given by Eq.|3]the width of the premelted boundary 
is predicted to be finite, with a value w m , at the melting point. 

In the MD results presented below we demonstrate the disordering behavior for two grain 
boundaries with relatively high and intermediate energy which we show can be quantitatively 



modeled by Eqs. [2] and [3j respectively. To compute quantitative values for the disjoining 
potentials from these simulations, we make use of the following relation between G(w) and 
the equilibrium distribution P(w) of grain-boundary widths sampled by the MD systems at 
a temperature T: P(w) oc exp [—AG(w)/kBT], where A is the interfacial area. Calculations 
of P(w) by MD, combined with an accurate knowledge of the bulk melting properties un- 
derlying AGf allows one to extract *f?(w), as described below. Since the MD simulations 
are based on a short-ranged classical interatomic-potential model, they do not include the 
long-ranged dispersion- force contributions to and thus sample ty sr (w) directly. 

METHODOLOGY 

In this study the simulations were based on an embedded-atom method (EAM) inter- 
atomic potential model for Ni developed by Foiles, Baskes and Daw (FBD) [39]. This po- 
tential was chosen as we have in previous work characterized both the bulk melting properties 
and solid-liquid interfacial free energies for FBD Ni; both sets of properties are prerequisites 
for a quantitative study of grain-boundary premelting. As reviewed in Ref. [20] the melting 
temperature for the FBD Ni potential has been bracketed to lie in the range T M =1709-1710 
K, and the value of the solid-liquid interfacial free energy has been calculated to be ^sl— 285 
mJ/m 2 [ID]. The latent heat for the potential was also calculated to be 0.015 eV/ A 3 . 

In the present study we considered the following three different grain boundary struc- 
tures: a E9 symmetric (011){411} 38.9° tilt boundary (hereafter referred to as the E9 tilt 
boundary), a Ell symmetric (011){311} 50.5° tilt boundary (hereafter referred to as the Ell 
tilt boundary), and a E9 (115) 120° twist boundary (hereafter referred to for convenience as 
the E9 twist boundary, even though we recognize that this boundary can also be described 
as a symmetric tilt boundary). These were selected from a large library of grain boundary 
structures that fit into moderately-sized periodic simulation cells, developed in the work of 
Olmsted et al.|41|. As shown in Table [IJ the selected boundaries have zero-temperature en- 
ergies (jgb) that span a range of values relative to 2jsl, and were thus expected to feature 
a range of premelting behavior. 
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Grain Boundary type 


L x (A) 


Ly (A) 


L z (A) 


Igb ijnJ/m 2 ) 




E9 (Oil) {411} 38.9° symmetric tilt 


237.86 


32.36 


31.68 


909 


1.5 


E9 (115) 120° twist 


253.23 


37.34 


38.80 


1440 


2.5 


Ell (011){311} 50.5° symmetric tilt 


246.55 


32.36 


33.02 


450 


0.79 



TABLE I: The grain boundary types used in the MD simulations along with the simulation cell 
sizes, reported at zero temperature. For reference the lattice constant of FBD Ni potential at 
zero temperature is 3.52 A. These boundaries span a large range of Jqb/^ISLi where Jq B is the 
zero-temperature grain-boundary energy, and jsl is the solid-liquid interfacial free energy. 
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FIG. 1: Computational cell used for molecular dynamics simulations at zero temperature. The 
boundary plane is in the y-z directions. The cell is periodic in y and z. 

Energy Minimization 



The optimization of grain-boundary structures at zero temperature made use of a simu- 
lation block with two grains separated by a flat grain boundary as shown in Fig [T] The cell 
was periodic in the plane of the boundary and non-periodic perpendicular to it. The grains 
were sandwiched between two slabs parallel to the grain-boundary plane. The atoms in each 
of these slabs were fixed relative to each other and could only undergo a rigid body motion. 
The slabs could move normal to the boundary, allowing volume expansion to maintain zero 
normal stress, and in the plane of the boundary, avoiding any resistance from the surfaces 
to translational movement of the grains relative to each other. Multiple trial configurations 
were built in order to minimize the boundary energy with respect to relative translations of 
the grains and with respect to the number of atoms in the grain boundary, as described in 
Ref. [41j. The energy of each trial configuration was minimized using the conjugate gradient 
method. After minimization, the energy of the grain boundary was computed as the total 
energy of the unconstrained atoms, less the bulk crystal energy for the same number of 
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atoms, divided by the area of the boundary. The dimensions of the computational cell used 
in these energy minimizations and in the subsequent MD simulations are given in Table [TJ 

Molecular-Dynamics Simulations 

Finite temperature simulations were performed by MD employing the LAMMPS code 
[42] . All simulations were performed with a constant number of atoms, constant volume 
and temperature (NVT ensemble). The temperature was maintained by using a Nose- 
Hoover thermostat [43] with a thermostat relaxation time of 0.1 ps, employing a time-step 
of 1 fs. The simulations began with the optimized zero-temperature geometry, which was 
equilibrated at different temperatures. Prior to each simulation at a given temperature, the 
periodic lengths parallel to the grain boundary plane were expanded according to the finite- 
temperature value of the lattice constant for the bulk crystal (determined separately). These 
periodic lengths were then held constant in all of the NVT MD simulations. The atoms in 
the slabs that were held fixed during the energy minimization procedure were made dynamic 
for the finite-temperature simulations, and the boundary conditions were modified such that 
both the grains terminated in free surfaces. 

As discussed in the next sub-section, two types of analysis were conducted on the simu- 
lated systems to study their premelting behavior. The first was a calculation of the excess 
volume as a function of temperature. In the simulations used for this analysis equilibration 
times were a few ns. The system volume was sampled over a total time of at least 10 ns at 
a frequency of 10 ps. 

The second analysis involved the calculation of equilibrium grain boundary width his- 
tograms. For these analyses, the simulations started at a given temperature with an equili- 
bration lasting 10 ns. Statistics were then obtained for the boundary width histograms from 
a total of 4000 snaphsots, sampled at a frequency of 10 ps. As discussed below, this number 
of snapshots and sampling rate ensured that at least a hundred statistically independent 
samples were obtained at each of the temperatures studied for the histogram analysis. 
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Calculation of Excess Volume and Width Histograms 

Prior to performing detailed analyses of width histograms, we performed an analysis of 
premelting behavior based on the temperature dependence of the excess volume. The excess 
volume was calculated by time averaging over values for single configurations, computed as 
follows. For a single snapshot the excess volume was computed from the distance between 
two lattice planes, one in each grain. These planes were chosen in each grain such that they 
were far from the boundary and the free surfaces. The volume of material that would lie 
between the two planes in a perfect crystal can be computed by counting the total number 
of atoms between the two planes (including half of the atoms in each of the two planes) and 
multiplying that number by the volume per atom of the bulk crystal at the same temperature 
(and zero pressure). The difference between the actual volume and the expected bulk volume, 
divided by the area of the boundary, is the excess volume. A slight linear dependence of 
this measured excess volume on the distance between the planes was found, presumably 
the result of the numerical error in the lattice constant at high temperature. However, any 
consistent choice was adequate for our purposes here, as the excess volume results were used 
mainly to determine the qualitative nature of the premelting behavior, as discussed below. 
The specific planes chosen for the excess volume were 1/4 and 3/4 of the way through the 
simulation cell. 

In order to compute equilibrium grain-boundary width distributions, P(w, T i ) ) for a given 
interface temperature T i) we proceeded as follows. For each snapshot the grain-boundary 
width w was determined using a scheme developed by Hoyt et al. [40] for the analysis 
of solid-liquid interface capillary fluctuations. In this approach each atom is assigned a 
structural order parameter, 4>i = j2 I ~ ^tj | 2 > w here are the actual positions of the 
12 nearest neighbors of atom % and r\- are the atom sites for the corresponding neighbors 
in the perfect crystal. The <^ values are then averaged in bins along the direction normal 
to the boundary and the point of inflection in the average order parameter profile is taken 
as the position of the grain boundary. In the case of grain boundaries two separate profiles 
are required. The first uses r\- for the crystal orientation of one of the two grains in the 
bicrystal, and for the second rf • is chosen based on the other grain. After these two order 



9 



parameter profiles (fa and fa) are obtained each is fit to the following function: 

fax) = a + 6tanh [(x — x )/d], (5) 

where a, 6, x and are fitting parameters, and the boundary width is derived from the 
difference in values of Xq determined from the fa and fa profiles. 

The analysis used to compute grain-boundary widths is illustrated graphically in Fig. [2] 
for a snapshot of the £9 twist boundary at an undercooling 2 K below the bulk melting 
temperature. In Fig. [2] the top panel shows the atoms, color coded according to the magni- 
tude of the difference in the order parameter values, lighter colors (red in the on-line figure) 
denoting atoms in an environment corresponding to either of the two grains, and darker 
colors (blue in the on-line figure) indicating an environment that does not correspond to 
either grain. The corresponding order parameter profiles with the fits to Eq. [5] (solid lines) 
are shown in the bottom of Fig. [2} The width is defined as the difference in the center 
positions of the fits to fa and fa. In this example, a width of 13.5 A is obtained by the 
analysis. 

It should be emphasized that the choice of an order parameter used to compute the 
grain boundary width is not unique (this point was discussed also in a recent related study 
by Williams and Mishin [33j). Other choices for the structural order parameter used to 
define an interface width would include the so-called centrosymmetry parameter [44J, the 
order-parameter used by Morris in studies of solid-liquid interfaces [45], excess mass used 
by Mellenthin[19j, a parameter introduced by Williams and Mishin [33j based on the local 
structure factor, and an order parameter introduced by von Alfthan et al. [34] based on 
structural units. Different choices for the order parameter are expected to lead to slightly 
different values for the interface width, and the quantitative values of the disjoining potential 
derived from them. For example, in the sharp-interface formalism we describe the limit of 
ty(w) for small w as the interfacial free energy of a hypothetical "dry" grain boundary - we 
would expect that this limiting value will depend on the way in which the width is defined, 
which is not unreasonable since with any choice of this definition a real boundary at zero 
temperature will have some small finite width. The ambiguities associated with the different 
choices for the definition of width is inherent in the use of a sharp-interface theoretical 
formalism to describe the properties of systems such as these with diffuse interfaces. In 
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FIG. 2: A snapshot of a premelted grain boundary at an undercooling of 2K and the corresponding 
average order parameter versus distance across the grain boundary. Atoms in the top panel are 
colored based on the order parameter value with lighter color (red in the on-line figure) representing 
atoms in either of the two grains and darker color (blue in the on-line figure) showing atoms in the 
grain boundary. The construction to determine the grain-boundary width w is illustrated in the 
bottom panel and a width of w « 13. 5A is obtained. 

applications of such sharp-interface models, however, the ambiguity presents no problem in 
practice provided that each of the bulk and interfacial thermodyamic quantities are defined 
consistently. 

In order to determine whether the simulation times employed in this study were sufficient 
to give adequate statistics for the width histograms, we estimated the correlation time for 
width fluctuations based on an analysis of the decay of a width autocorrelation function 
C(At) defined as follows: 



C(At) = (w(t)w(t + At)) - {wf 



(6) 



In Eq. [6j w(t) denotes the instantaneous width of the grain boundary at a time £, and (...) 
denote ensemble (time) averages. The time required for the autocorrelation function to 
decay by a fraction of 1/e was taken as a measure of the correlation time. 

The correlation times were found to vary significantly with boundary type and temper- 
ature. The maximum value of the correlation time was obtained as 100 ps for the S9 tilt 
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boundary at a temperature 2 K below Tm- Given that 4000 snapshots were sampled in 
the MD simulations, at a frequency of 10 ps, at least a hundred independent samples were 
obtained for each of the histograms presented below. 

Calculation of Disjoining Potentials from Width Histograms 

As discussed in the previous section, the probability of observing a given boundary width 
at an interface temperature is related to the total free energy of the system as follows: 

P(w,T z ) = dexpi-AGiw^/knTii (7) 

where G{w ) T i ) is given by Eq. [I] In the above expression the subscript i denotes one of the 
seven different temperatures for which the width histograms have been calculated by MD: 
Ti = 1710, T 2 = 1708, T 3 = 1705, T 4 = 1700, T 5 = 1690, T 6 = 1680, T 7 = 1650X, which 
represent undercoolings relative to the bulk melting temperature, T M =1710 K, ranging 
between zero and 60 K. Equation [7] emphasizes the fact that each undercooling yields a 
histogram that spans a different range of widths, and the unknown normalization constants 
at each temperature are denoted as C{. 

As discussed in the next section, analyses of the calculated width histograms were un- 
dertaken to compute disjoining potentials for two different boundaries, one of which (the 
£9 twist) featured a diverging interface width as T M is approach from below, and another 
(the S9 tilt) whose width remains finite at T M - In the former case the disjoining potential 
was modeled using the form given by Eq. [2j while for the latter use was made of Eq. [3| 
The unknown parameters in these expressions can be obtained in two ways based on the 
MD-calculated width histograms using Eq. [7| The first is to fit the histogram data at each 
temperature independently. This gives a range of values for the disjoining potential param- 
eters, which are expected to be more accurate at the higher temperatures where one has 
access to the largest range of width values. A refined estimate of the potential parameters 
can be obtained from a histogram analysis using all of the data at once, as follows. 

From Eq. [I] the disjoining potential can be written as 

*(w) = G(w,T l )-wAG f (T l ) (8) 
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where AGf is the bulk free energy difference between liquid and solid per unit volume. For 
the EAM Ni system studied here the temperature dependence of AGf is known accurately 
from previous studies. For the small range of temperatures considered in the MD simulations 
of width histograms, which are close to the bulk melting temperature, AGf can be accurately 
computed from the following expression: 

AG f = L(ATi/T M ) (9) 

where L is the latent heat per unit volume, and ATi = T M — Ti is the interface undercool- 
ing. From Eq. [7| the disjoining potential can be written in terms of the equilibrium width 
distribution as: 

W(w) = -{k B T i )lnP{w ) T i )/A i - AGfW + ai (10) 

where the are unknown constants related to C$ in Eq.[7j These parameters lead to constant 
offsets when the quantity -fc^T^ In P(w, Tj) jAi — AGfW is plotted versus w for each interface 
temperature. To construct the disjoining potential a least squares fit is used to refine the 
values of the shift parameters (a$) along with the potential parameters (A7 and S in Eq. [2] 
or Ai, 61, A 2 , and S 2 in Eq. [3]). In this procedure the initial values for the shift parameters 
were estimated "by eye" to give maximal overlap between the data, and the initial values 
for the potential parameters were obtained from an average of the independent fits to the 
data at separate temperatures (described above). The values of the parameters were then 
iteratively refined to minimize the square of the differences between the MD data and the 
analytical expressions for ty(w). 

Effects of Thermostat, System Size and Boundary Conditions 

In the course of this work additional studies were also conducted to investigate the effects 
of the details of the simulation methodology on the resulting disjoining potentials. These 
additional analyses were performed for the £9 twist boundary in Ni as well as a series of 
(100) tilt boundaries in Ni and bcc Fe. 

The first issue addressed in these additional simulations was the effect of thermostat. 
Equation [7] assumes that the measured grain-boundary widths sample a canonical ensem- 
ble, and in order to achieve this by MD simulations we employed a standard Nose-Hoover 
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thermostat [43j in the present work. An alternative choice would be to use a Langevin 
thermostat [46j , which would relax temperature fluctuations faster locally and produce dif- 
ferent dynamics in the system, but would be intended to sample the same ensemble. To 
test for unexpected effects of dynamics on the equilibrium width histograms and disjoining 
potential, two boundaries showing continuous premelting behavior (i.e., diverging widths as 
the melting temperature is approached from below), one in bcc Fe and one in fee Ni, were 
simulated using both a Nose-Hoover thermostat and a Langevin thermostat with a time 
constant of 0.25 ps. No statistically significant differences in the disjoining potential were 
found with the two thermostats. Further, we tested a range of reasonable time constants 
for both thermostats, finding that statistically significant effects on the calculated disjoining 
potentials were again absent unless the relaxation time for the Langevin thermostat was 
reduced to extremely short time scales, on the order of the inverse of the vibrational fre- 
quencies in the system. Interestingly, in simulations which used the same time constant of 
0.25 ps, the correlation time for width fluctuations was found to be substantially reduced in 
the simulations with the Langevin relative to those with the Nose-Hoover thermostat. The 
Langevin thermostat thus offers the potential advantage of providing improved statistics 
related to width histograms for a given simulation time. 

The second issue investigated was the effect of system size on the calculated disjoining 
potential. For atomically rough grain boundaries, such as those investigated in this work, 
previous theoretical studies have demonstrated that capillary fluctuations should give rise 
to a weak dependence of the disjoining potential on cross-sectional area [IU [28]. To check 
for any unexpected large system-size effects, simulations were performed for the S9 twist 
boundary using both the cross-sectional area given in Table [TJ as well as a value that was 
four times larger (i.e, with periodic lengths that were doubled in each of directions parallel 
to the boundary). Although the width histograms obtained for the larger system were much 
narrower, as expected from Eq. [7| it was still possible to obtain data over a sufficiently large 
range of widths to extract a disjoining potential. The results for the disjoining potential were 
found to be consistent with those obtained from the smaller system provided that a small 
effect of system size on the bulk melting temperature (amounting to a roughly one-degree 
lowering of T M for the larger system, which is within the accuracy of the known value of 
the bulk melting point) was accounted for in the analysis. Similar results were obtained in 
analyses of system size effects on calculated disjoining potential for a high-energy (100) tilt 
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boundary in Ni. 

A final issue that was investigated in these studies concerns the effect of boundary condi- 
tions on the calculated width distributions. For a high-angle (100) symmetric tilt boundary 
in Fe this was investigated by computing width distributions at 10 K undercooling using 
both the free-surface boundary conditions employed in the current study, as well as a full pe- 
riodic boundary condition in a simulation cell containing two grain boundaries and a periodic 
length normal to the boundaries that was allowed to be dynamic to maintain zero normal 
stress. No statistically significant differences were found between the width distributions 
derived with the two different boundary conditions. 

RESULTS 

The effect of temperature on the structure of the three boundaries considered in this 
work is illustrated in Figs. [3]and[4| Figure [3] plots the calculated excess volume versus the 
logarithm of the interface undercooling, and displays three qualitatively different types of 
behavior. The £9 twist boundary has an excess volume (represented by the red circles) 

T 1 1 I I I I | 1 1 1 1 I I I I I 1 1 1 1 I I I I | 

• 19 <1 15> 120° twist 

♦ £9 symmetric <01 1> {41 1 } 38.9°tilt 
a £1 1 symmetric <01 1> {31 1 } 50.5° tilt _ 

• 

* ♦ ♦ ♦ 

♦ • 

* • • 

A AAA . 

A 1 * t " 

J I I I I I I I I I I 

10 100 1000 

Undercooling (K) 

FIG. 3: Excess volume versus temperature plot for three different grain boundaries in Nickel. The 
red circles correspond to a E9 (115) twist grain boundary, the green diamonds correspond to a £9 
symmetric (011) {411} tilt grain boundary and the black triangles correspond to a £11 symmetric 
(011) {311} tilt grain boundary. The three boundaries show very different behavior. 

that continuously increases as the melting temperature is approached from below. The 
behavior displayed by this boundary in Fig. [3] is consistent with a logarithmic divergence of 
the interface width, as would be expected if the disjoining potential is of the form given by 
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Eq. [2j At the opposite extreme, the Ell tilt boundary (represented by the black triangles) 
has an excess volume that is relatively small and depends only weakly on temperature. The 
E9 tilt boundary (represented by the green diamonds) displays an intermediate behavior. 
The excess volume rises with increasing temperature at a rate that initially tracks that of 
the continously premelting E9 twist boundary. At high temperature, the rate of increase of 
the excess volume decreases and the boundary maintains a finite excess volume at T M . 

Representative snapshots and calculated width histograms are shown in Fig. [| for each of 
the three boundaries at 1708 K (two degrees below T M )- The width histogram corresponding 
to the E9 twist boundary is much broader than those for the other two boundaries. This 
boundary samples relatively large widths, and P(w) shows pronounced asymmetry with 
an extended tail to large w values. The Ell tilt boundary, by contrast, features a width 
histogram that is very narrow and is centered on relatively small values of w. The width 
distribution for this boundary is roughly Gaussian in shape, with no detectable asymmetry 
towards large values of w. The width distribution for the E9 tilt boundary shows features 
intermediate between the other two. The average width is larger than that of the Ell 
boundary, and the width distribution is considerably broader. Compared to the E9 twist 
boundary, however, the asymmetry and the tail extending to larger widths is not nearly as 
pronounced in the width distribution for the E9 tilt boundary. 

The MD snapshots in Fig. [4] next to each histogram provide additional insights into the 
nature of the high-temperature structural disorder present in each of the three boundaries. 
As in Fig. [2| the atoms in these snapshots have been colored based on the difference between 
values of the structural order parameters (0i and 2 ) defined above. Lighter colors (red in the 
on-line figure) indicate atoms with an environment corresponding to one of the two grains, 
while darker color (blue in the on-line figure) denotes a disordered environment distinct 
from those in either of the grains. The snapshots and width histograms show that the E9 
twist boundary at 1708 K displays thick premelted layers, consistent with the continuous 
premelting behavior suggested by the excess volume. The Ell tilt boundary is seen to 
display appreciable disorder only over a width that is on the scale of one atomic plane; 
this boundary is observed in the MD simulations to remain highly ordered up to Tm- The 
intermediate behavior of the E9 tilt boundary is characterized by disorder that extends over 
a distance of several atomic planes. An important feature of the disorder observed in this 
boundary is the nonuniform character of the widths along the area of the boundary. As 
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FIG. 4: The distribution function P{w) vs w from the MD simulations for the three different 
boundaries at 1708 K. The green diamonds and the red circles correspond to the £9 tilt and twist 
grain boundaries respectively. TheSll tilt boundary is represented by the black triangles. The 
lines represent a least square fit of the premelting models of Eqs. [I] and [2] for the £9 twist grain 
boundary, Eqs. [I] and [3] for the £9 tilt boundary and a gaussian fit for the £11 tilt boundary. 
The snapshots next to each figure correspond to the widest grain boundary for each boundary and 
are color coded according to the fa values. Lighter color (red in the on-line figure) indicates an 
atom with an environment corresponding to one of the two grains, darker color (blue in the on-line 
figure) represents atoms in some other environment. 

illustrated in the snapshot of this boundary in Fig. [4j thick and extended regions of disorder 
up to a nanometer or more in thickness are often observed in part of the boundary, with 
much narrower, more ordered regions in between. For the remainder of this section we will 
focus on the behavior of the two £9 boundaries, showing that the temperature dependence 
of the width histograms can be accurately described by the disjoining-potential formalism 
described above. 




FIG. 5: Snapshots from a MD simulation at an undercooling 2 K illustrating the dynamic nature of 
the grain-boundary width. The snapshots corresponding to the £9 tilt grain boundary are on the 
right of the panel. The snapshots to the left are from the £9 twist boundary discussed elsewhere 
[20] . The atoms are colored by the difference between the two order parameters as described above. 
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FIG. 6: The distribution function P(w) vs w from the MD simulations. The function corresponding 
to the £9 tilt grain boundary is in the right of the panel. The function to the left is from the £9 
twist boundary discussed elsewhere [20] . The lines represent least square fits of the premelting 
model of Eqs. [I] and [2] for the £9 twist grain boundary and Eqs. [I] and [3] for the D9 tilt boundary. 

Figure [5] further illustrates the nature of the disorder present in the S9 twist and tilt 
boundaries. The snapshots were obtained from a 40 ns MD simulation at 1708 K, and 
represent the largest, smallest and an average value for the interface widths. Results for the 
£9 twist boundary are shown on the left and those for the £9 tilt boundary are on the right 
of Fig. [5j The snapshots illustrate the highly dynamic nature of the boundary structures at 
1708 K. The large fluctuations in interface width displayed by these boundaries forms the 
basis for the histogram analysis underlying calculations of the disjoining potential. 

The calculated width histograms for both £9 boundaries are displayed over a range of 
interface undercoolings in Fig. [6} The histograms on the left and right in Fig. [6] correspond 
to the £9 twist and tilt boundaries, respectively. For both of the boundaries the histograms 
become broader as the bulk melting temperature is approached (i.e., with decreasing un- 
dercooling). For the £9 twist boundary the system is observed to melt completely over the 
course of the 20 ns runs at the bulk melting temperature, and histograms can be obtained 
for this boundary only for finite values of the undercooling. By contrast, the width of the 
£9 tilt boundary remains finite at T M and the width histogram can be calculated at zero 
undercooling for this boundary, as shown in Fig. [HJ The inset in each of the panels shows 
the quantity — fcgT^log [P(w,Ti)]/A — AGfW versus w for all of the interface temperatures. 
The data for the £9 tilt boundary shows clear evidence of an attractive minimum in the 
disjoining potential, as will be discussed further below. 
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The solid lines in Fig. [^represent the least square fits of the MD data for P(w) to the 
disjoining-potential formulas given by Eqs. [7| [T] and [2] for the £9 twist boundary, and Eqs.[7| 
[T] and [3] for the £9 tilt boundary. For the £9 twist boundary the potential parameters 
obtained from the separate fits to the data for the individual temperatures span the range 
5=0.25 to 0.29 nm, and A7=101 to 150 mJ/m 2 . For the £9 tilt boundary the fitted values 
for Si and 82 ranged between 0.142 to 0.144 and 0.143 to 0.144 nm, respectively, Ai — A 2 
spanned the values 103 to 163 mJ/m 2 , and A 2 /Ai took values between 1.003 and 1.006. 
With these fitted potential parameters, the disjoining potential for the £9 tilt boundary 
exhibits a weak minimum at a width w m with values ranging between 0.32 and 0.36 nm, 
and a depth relative to 2jsl varying between -7 and -11 mJ/m 2 . The fact that the MD 
calculated width histograms can be accurately described by disjoining potentials with a 
relatively narrow range of fitted potential parameter values indicates that the formalism 
described in the previous section represents a valid model for describing the temperature 
dependence of the structural disorder observed in these boundaries. 

In order to refine the calculation of the disjoining potential for the £9 twist and tilt 
boundaries, we employ the histogram analysis described in the previous section, involving 
a refinement of both the shift parameters and the potential parameters in Eqs. [2] and [3] 
for the twist and tilt boundaries, respectively. The resulting fits are shown in Fig. [7] and 
correspond to the following values for the potential parameters: S = 0.25 nm and A7 = 156 
mJ/m 2 for the £9 twist boundary, and A 1 - A 2 = 103 mJ/m 2 , A 2 /Ai = 1.003, Si = 0.1471 
nm, S 2 = 0.1474 nm for the £9 tilt boundary. These parameter values are consistent with 
the values given above from the independent fits. The excellent agreement of the fits with 
the MD data in Fig. [7] again indicates that the disjoining-potential formalism represents 
an accurate framework for modeling the premelting behavior of these £9 boundaries. The 
analysis used to obtain the results in fits in Fig. [6] assumed a melting point of T M — 1710 
K. If the melting temperature is changed even by one degree, i.e., T M = 1709X then poor 
fits to P(w) are obtained for the data at the lowest undercoolings. 

The calculated disjoining potentials in Fig. [7] are characterized by the following features. 
For the £9 tilt boundary the disjoining potential has a minimum at a finite width, w m , which 
corresponds to the average equilibrium grain boundary width at the melting temperature. 
The potential is repulsive for w < w m and attractive for w > w m . In contrast, the £9 twist 
boundary features a purely repulsive disjoining potential that is well modeled by exponential 
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FIG. 7: An illustration of the histogram method used to extract the disjoining potential. It shows 
the merged data from individual histogram data used to reproduce the complete disjoining potential 
ty(w) for the £9 twist and the D9 tilt boundaries. The solid line is the best fit to the exponential 
decay given in Eq. [2] for the S9 twist boundary, and a double exponential function in Eq. [3] for the 
£9 tilt grain boundary. 

form given in Eq. [2j as expected based on the logarithmic divergence of the excess volume 
calculated for this boundary (see Fig. [3]). 



SUMMARY AND DISCUSSION 



In this paper we have presented a detailed description of a method for computing the 
disjoining potential for grain-boundary premelting and grain coalescence from an analysis of 
width fluctuations measured in equilibrium MD simulations. The approach has been applied 
to two grain boundaries in an EAM model of elemental Ni. For the £9 twist boundary, 
which features an excess volume that increases logarithmically as T M is approached from 
below, the measured width histograms are consistent with a disjoining potential that decays 
exponentially, with a decay length of approximately 0.25 nm and a maximum value at zero 
width of approximately 160 mJ/m 2 relative to 2^sl- For the S9 tilt boundary, which features 
an excess volume that remains finite at T M) the measured width histograms are shown to be 
consistent with a disjoining potential that features a weak attractive minimum at w m ^0.34 
nm, with a depth relative to of approximately -8 mJ/m 2 . 

It is interesting to compare the results of the present work with previous studies of grain- 
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boundary premelting based on MD simulations, phase-field theory and PFC calculations. 
The exponential form of the disjoining potential for the £9 twist boundary presented here 
and in Ref. [20j is qualitatively consistent with the results of previous investigations where 
a diverging grain-boundary width has been observed for high-energy boundaries in a variety 
of different systems (e.g., [25l [26l [28j [34)). The disjoining potential calculated here for the 
S9 tilt boundary features a weak attractive minimum at w m , and is replusive and attractive 
for smaller and larger values of w, respectively. This form for the disjoining potential 
corresponds to a grain-boundary whose width increases with T at low temperature, but 
remains finite at and above T M (up to some maximum superheating) . This type of behavior 
for the temperature dependence of the grain boundary width is qualitatively similar to that 
found in recent MD studies for one of three twist boundaries in Si [34] and a symmetric tilt 
boundary in Cu [33]. 

As discussed in the introduction section disjoining potentials with attractive minima, 
qualitatively similar to that calculated here for the £9 tilt boundary, were obtained in PFC 
calculations by Mellenthin et al.[19] for low-angle tilt boundaries in a model two-dimensional 
hexagonal system. The boundaries which displayed this type of disjoining potential had 
highly non-uniform interface structures, containing premelted regions localized on grain- 
boundary dislocation cores and separated by well ordered regions which we will refer to 
as solid "bridges". For these structures, the width corresponding to the minimum in the 
disjoining potential, which represents an average over the area of the boundary, corresponds 
to the thickness of an equivalent uniform layer with the same amount of "liquid" as contained 
in the premelted cores. This area-averaged width can be quite small (e.g., on the order of 
the atomic spacing) even though the radius of the premelted cores is much larger. 

It is interesting to compare the structures observed in the PFC calculations of Mellenthin 
et al. [19J with those for the £9 tilt boundary considered in the present study. Visual 
inspection of the numerous snapshots showed that the disorder in this boundary is indeed 
often highly non-uniform, as illustrated in Fig. |8j 

This figure shows a representative snapshot from a simulation at T = 1708 K, viewed 
down the tilt axis. In Fig. |8j the structural disorder in the boundary is clearly nonuniform 
- the region of high structural disorder is highlighted by the grey ellipse, while the red lines 
connected across the boundary plane highlight the ordered solid bridges. Thus, even though 
the £9 boundary studied in the present work has a misorientation angle (38.9 degrees) 
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FIG. 8: Snapshot of the £9 tilt boundary at 1708 showing the non-uniform structural disorder 
in the grain boundary. The grey ellipse shows the regions with high disorder while the red lines 
highlight the ordered bridges. 

which is far too large for its structure to be described in terms of separated grain-boundary 
dislocation cores, the non-uniform nature of the structural disorder and the presence of solid 
"bridges" is qualitatively very similar to the types of structures observed in Ref . [19] . 

It should be emphasized, however, that the structure of the £9 boundary observed in 
the MD simulations is highly dynamic (c.f., Fig. [5]) such that the solid bridges form and 
disappear rapidly on the time scale of the simulations. This behavior could have interesting 
consequences for the shear response of such boundaries, as the solid bridges are expected to 
offer enhanced resistance to shear which would otherwise be expected to be very limited for 
a premelted grain boundary (e.g., [26]). The shear response of such boundaries would thus 
be an interesting topic for future MD studies. 

ACKNOWLEDGMENTS 

Work at Northeastern, UC Davis and Berkeley was supported by the Director, Office of 
Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division, of 
the U. S. Department of Energy (DOE), under contracts DE-FG02-07ER46400, DE-FG02- 
06ER46282, and DE-AC02-05CH11231. JJH acknowledges financial support from a Natural 
Sciences and Engineering Research Council (NSERC) of Canada Discovery grant. All the 



22 



authors acknowledge support from the DOE Computational Materials Science Network pro- 
gram. This research used resources of the National Energy Research Scientific Computing 
Center, which is supported by the Office of Science of the U.S. Department of Energy under 
Contract No. DE-AC02-05CH11231. We also acknowledge helpful discussions with R. G. 
Hoagland. 



[1] M. Rappaz, A. Jacot, and W. J. Boettinger, Met. Mater. Trans. 34A, 467 (2003) 

[2] M. Rappaz, J.M. Drezet, and M. Gremaud, Met. Mater. Trans. 30A, 467 (1999) 

[3] N. Wang, S. Mokadem, M. Rappaz and W. Kurz, Acta Mater. 52, 3173 (2004) 

[4] M. Asta, C. Beckermann, A. Karma, W. Kurz, R. Napolitano, M. Plapp, G. Purdy, M. Rappaz 

and R. Trivedi, Acta Mater. 57, 941 (2009) 
[5] J. G. Dash, H. Fu and J.-S. Wettlaufer, Rep. Prog. Phys. 58, 115 (1995) 
[6] R. W Balluffi and R. Maurer, Scripta Metall. 22, 709 (1988) 
[7] R. A. Masamura, M. E. Glicksman and C. L Void, Scripta Metall. 6, 943 (1972) 
[8] T. Watanabe, S. I. Kimura, S. Karashima, Philos. Mag. A. 49, 845 (1984) 
[9] C. L. Void and M. E. Glickman, in: The Nature and Behavior of Grain Boundaries, edited 

by H. Hu, Metallurgical Society of AIME, Plenum Press, New York, 171-183 (1972) 
[10] F. Inoko, T. Muraga, T. Nakano, Y. Yoshikawa, T. Interface Sci. 4, 263 (1997) 
[11] S. V. Divinski, M. Lohman, C. Herzig, C Straumal, B. Baretzky, W. Gust, Phys. Rev. B 71, 

104104 (2005) 

[12] J. Lou, V. K. Gupta, D. H. Yoon and H. M. Meyer, Applied Physics Letters 87, 231902 (2005) 
[13] V. K. Gupta, D. H. Yoon and H. M. Meyer, J. Lou, Acta. Mater. 55, 3131 (2007) 
[14] A. E. Lobkovsky and J. A. Warren, Physica D 164, 202 (2002) 

[15] M. Tang, W. C. Carter and R. M. Cannon, Phys. Rev. Lett. 97, 075502 (2006); Phys. Rev. B 
73, 24102 (2006) 

[16] Y. Mishin, W. J. Boettinger, J. A. Warren and G. B. McFadden, Acta Mater. 57, 3771 (2009) 
[17] N. Wang, R. Spatschek, and A. Karma, "Multi-phase- field analysis of short-range forces be- 
tween diffuse interfaces", arXiv:0912.4627vl [cond-mat.mtrl-sci] (2009) 
[18] J. Berry, K. R. Elder and M. Grant, Phys. Rev. B, 77, 224114 (2008) 
[19] J. Mellenthin, A. Karma and M. Plapp, Phys. Rev. B, 78, 184110 (2008) 



23 



[20] J. J. Hoyt, D. Olmsted, S. Jindal, M. Asta and A. Karma, Phys. Rev. E 79, 020601 (R) (2009). 

[21] T. Nguyen, P. S. Ho, T. Kwok, C. Nitta, and S. Yip, Phys. Rev. Lett. 57, 1919 (1986) 

[22] P. S. Ho, T. Kwok, T. Nguyen, C. Nitta, and S. Yip, Scr. Metall. 19, 993 (1985) 

[23] G. Ciccotti, M. Guillope, and V. Pontikis, Phys. Rev. B 27, 5576 (1983) 

[24] M. Guillope, G. Ciccotti, and V. Pontikis, Surf. Sci 144, 67 (1984) 

[25] J. Q. Broughton and G. H. Gilmer, Phys. Rev. Lett. 56, 2692 (1986) 

[26] J. Q. Broughton and G. H. Gilmer, Modelling. Simul. Mater, sci. Eng. 6, 87-97 (1998) 

[27] J. F. Lutsko, D. Wolf, S. Yip, S. R. Phillpot and T. Nguyen, Phys. Rev. B 40, 2841 (1989) 

[28] R. Kikuchi and J. W. Cahn, Phys. Rev. B, 21, 1893 (1980) 

[29] A. Suzuki, Y. Mishin, Mater. Sci. Forum 502, 157 (2005) 

[30] G. Besold, O. G. Mouritsen, Comp. Mater. Sci. 18, 225 (2000) 

[31] P. Keblinski, D. Wolf, S. R. Phillpot, H. Gleiter, Philos Mag A 79, 2735 (1999) 

[32] J. Lu and J. A. Szpunar, Interface Sci. 3, 143 (1995) 

[33] P. L. Williams, Y. Mishin, Acta. Mater. 57, 3786 (2009) 

[34] S. Von Alfthan, K. Kaski, A. P. Sutton, Phys. Rev. B 76, 245317 (2007) 

[35] R. Lipowsky, Phys. Rev. Letter 57, 2876 (1986) 

[36] D. R. Clarke, J. Am. Ceram. Soc, 70, 15 (1987) 

[37] B. Widom, J. Chem. Phys., 68, 3878 (1978) 

[38] R. Lipowsky and M. E. Fisher, Phys. Rev. B 36, 2126 (1987) 

[39] S. M. Foiles, M. I. Baskes and M. S. Daw, Phys. Rev. B, 33, 7983 (1986) 

[40] J. J. Hoyt, M. Asta and A. Karma, Mater. Sci. Engr.,R41, 121 (2003) 

[41] D. L. Olmsted, Acta Mater, 57, 2793 (2009) 

[42] S. J. Plimpton, J. Comp. Phys. 117, 1-19 (1995) 

[43] S. Melchionna, G. Ciccoti, and B. L. Holian, Mol. Phys. 78, 533 (1993) 

[44] C. L. Kelchner, S. J. Plimpton, J. C. Hamilton, Phys. Rev. B 58 11085 (1998) 

[45] J. R. Morris, Phys. Rev. B 66, 144104 (2002) 

[46] Dunweg and Paul, Int J of Modern Physics C, 2, 817-27 (1991) 



24 



